#ifndef GADGET_SNAPSHOT_IO_H
#define GADGET_SNAPSHOT_IO_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "utils/utils.02.h"
#include "Timer/timer.01.h"
#include "FileNameDecomposer/fnd.02.h"
#include "GadgetParticleType/gpt.01.h"


class GSIO {
	public:
	
	FND fnd;	
		
	struct io_header_1 {
		
		int      npart[6];
		double   mass[6];
		double   time;
		double   redshift;
		int      flag_sfr;
		int      flag_feedback;
		int      npartTotal[6];
		int      flag_cooling;
		int      num_files;
		double   BoxSize;
		double   Omega0;
		double   OmegaLambda;
		double   HubbleParam; 
		char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  /* fills to 256 Bytes */
	
	} header1;
	
	int NumPart, Ngas;
	
	particle_data *P, *_p[6];
	
	uint *Id;
	uint blocksize;
	uint pTypeNum;
	char *mem;
	char *filename, *filename_base;
	uint *ofs;
	
	class TN {
		public:
			int N[7];
			int &operator =(int i)   {  N[0]=i; return i;  }
			int &operator [](int i)  {  return N[i+1];     }
			operator int(){	return N[0]; }
	} N;
	
	GSIO():
	P(0),
	Id(0),
	blocksize(4),
	pTypeNum(6),
	mem(0),
	filename(0),
	filename_base(0),
	ofs(0)
	{ }
	
	~GSIO(){
		if (P) delete [] P;
		if (mem) delete [] mem;
		if (Id) delete [] Id;
	}

	double Time, Redshift;
	double readtime;
	double readspeed;
	uint filesize;
	int N_withVarMassess;
	
	char* readraw(const char* filename) {
		
		filesize = getFileSize(filename);
		char *mem = new char[filesize];
		
		Timer tr; 
		FILE *fs = fopen(filename, "r");
		if (fs==NULL) {printf("\033[22;31merror: can't read the file\033[22;35m '%s'\033[22;30m\n", filename); return NULL;}
		
		tr.set();
		fread(mem, 1, filesize, fs); 
		tr.get();
		
		readspeed = filesize/tr.sin;
		
		fclose(fs);
		
		return mem;
		
	}

	char* getSplitFileName(const char *basefilename, uint ifile){
		char *splitfilename = new char[256];
		sprintf(splitfilename, "%s.%u", basefilename, ifile);
		return splitfilename;
	}
	
	void read(const char* filename){
		
		Timer tr; 
		
		fnd.decompose(filename);
		
		if (this->filename) delete [] this->filename;
		this->filename = new char[256];
		sprintf(this->filename, "%s", filename);
		
		uint ifile=0;
		
		char *mem0;
		
		uint *ofs = new uint[pTypeNum];
		memset(ofs,0,pTypeNum*sizeof(uint));
		
		uint filesize(0);
		
		do {
			
		
		if (mem)	delete [] mem;
		tr.set();
		mem0 = mem = readraw(this->filename);
		
		tr.get();
		filesize+=this->filesize;
		
		//header
		mem+=blocksize;
		memcpy(&header1, mem, sizeof(io_header_1)); 
		mem+=sizeof(io_header_1);
		mem+=blocksize;
		
		
		// local
		uint NumPart=0;
		uint Ngas= header1.npart[0];
		for (uint i=1; i<pTypeNum; i++)
		NumPart+= header1.npart[i];
		uint N=NumPart+Ngas;
		
		// global
		this->NumPart=0;
		this->Ngas= header1.npartTotal[0];
		for (uint i=1; i<pTypeNum; i++)
		this->NumPart+= header1.npartTotal[i];
		
		
		this->N=this->NumPart+this->Ngas;
		for (uint i=0; i<pTypeNum; i++)
		this->N[i] = header1.npartTotal[i];
		
		
		
		// offsets
		if (ifile==0) {
			if (this->ofs) delete [] this->ofs;
			this->ofs = new uint[pTypeNum];
			memset(this->ofs,0,pTypeNum*sizeof(uint));
			for (uint k=0; k<pTypeNum-1; k++)
			this->ofs[k+1] = this->ofs[k] + header1.npartTotal[k];
		}

		// P
		if (!ifile) {
			if (P) delete [] P;
			P = new particle_data[this->N];
		}

		_p[0] = P;
		for (uint i=1; i<pTypeNum; i++)
		_p[i] = _p[i-1] + this->N[i-1];



		// Pos
		mem+=blocksize;
		for (uint k=0; k<pTypeNum; k++)
		for (uint j=0; j<header1.npart[k];j++) {
		memcpy(P[ this->ofs[k]+ofs[k]+j ].Pos, mem, 3*sizeof(float));  
		mem+=3*sizeof(float); 
		}
		mem+=blocksize;
														
		// Vel
		mem+=blocksize;
		for (uint k=0; k<pTypeNum; k++)
		for (uint j=0; j<header1.npart[k];j++) {
		memcpy(P[ this->ofs[k]+ofs[k]+j ].Vel, mem, 3*sizeof(float));  
		mem+=3*sizeof(float); }
		mem+=blocksize;
		
		// Id
		if (!ifile) {
			if (Id) delete [] Id;
			Id = new uint[this->N];
		}
		
		// Ids
		mem+=blocksize;
		for (uint k=0; k<pTypeNum; k++) {
		memcpy(&Id[this->ofs[k]+ofs[k]], mem, header1.npart[k]*sizeof(uint));  
		mem+=header1.npart[k]*sizeof(uint);
		}
		mem+=blocksize;

		// Mass
		N_withVarMassess=0; // lets check whether the var masses block is present
		for (uint k=0; k<pTypeNum; k++)
		if (header1.mass[k]==0)
		N_withVarMassess+=header1.npart[k];
		
		
		
		if (N_withVarMassess>0) mem+=blocksize; // do it if the var masses block is present
		for (uint k=0; k<pTypeNum; k++) 
		for (uint i=0; i<header1.npart[k]; i++) {
			if (header1.mass[k]==0) { // we read the block
				memcpy(&P[this->ofs[k]+ofs[k]+i].Mass, mem, 1*sizeof(float));  
				mem+=1*sizeof(float); 
			}
			else // in this case we use common masses
				P[this->ofs[k]+ofs[k]+i].Mass = header1.mass[k];
			
		}
		if (N_withVarMassess>0) mem+=blocksize;

		
		// Gas
		if(header1.npart[0]>0) {
			
			// Internal Energy
			mem+=blocksize;
			for (uint j=0; j<header1.npart[0]; j++) {
			memcpy(&P[*ofs+j].U, mem, 1*sizeof(float));  
			mem+=1*sizeof(float);	}
			mem+=blocksize;

			// Density
			mem+=blocksize;
			for (uint j=0; j<header1.npart[0]; j++) {
			memcpy(&P[*ofs+j].Rho, mem, 1*sizeof(float));  
			mem+=1*sizeof(float);  }
			mem+=blocksize;
		
			
			if(header1.flag_cooling) { // this is for the cooling in the user guid it says it is currently unused
				mem+=blocksize;
				for (uint j=0; j<header1.npart[0]; j++) {
					memcpy(&P[*ofs+j].Ne, mem, sizeof(float));
					mem+=sizeof(float);
				}
				mem+=blocksize;
			} else 
				for (uint j=0; j<header1.npart[0]; j++)
				P[*ofs+j].Ne= 1.0;
		
		}
		
		
		for (uint k=0; k<pTypeNum; k++) 
		ofs[k]+=header1.npart[k];
		
		mem = mem0;	
		
		if (header1.num_files>1) { 
			
			ifile++; 
			
			if (this->filename) delete [] this->filename;
			this->filename = getSplitFileName(fnd.fileName.c_str(), ifile);
			
		} else break;
		
		} while (ifile<header1.num_files);
 
		Time= header1.time; // read snapshot time (or exspansion factor) 
		Redshift= header1.time; // read the snapshot redshift 
		
		readspeed=filesize/tr.sin_tot;
		readtime=tr.sin_tot;
		this->filesize=filesize;
		delete [] ofs;


	}
	
	void write(const char *filename){
		
		char *mem0(0);
		
		uint *ofs = new uint[pTypeNum];
		memset(ofs,0,pTypeNum*sizeof(uint));
		
		// I need offset
		if (this->ofs) delete [] this->ofs;
		this->ofs = new uint[pTypeNum];
		*(this->ofs)=0;
		for (uint k=0; k<pTypeNum-1; k++)
		this->ofs[k+1] = this->ofs[k] + header1.npartTotal[k];
		
		for (uint ifile=0; ifile<header1.num_files; ifile++){
		
			for (uint k=0; k<pTypeNum; k++) 
			if (header1.npart[k]) {
				header1.npart[k]=header1.npartTotal[k]/header1.num_files;
				if (ifile==header1.num_files-1) header1.npart[k]+=header1.npartTotal[k]%header1.num_files;
			} else header1.npart[k]=0;
		
			
			uint n_tot(0);
			for (uint k=0; k<pTypeNum; k++) 
			if (header1.npart[k]) n_tot+=header1.npart[k];
		
			uint n_totvm(0);
			for (uint k=0; k<pTypeNum; k++) 
			if (header1.npart[k] && !header1.mass[k]) n_totvm+=header1.npart[k];
			
			//sizeoffile 
			filesize = sizeof(io_header_1) +
			n_tot*sizeof(particle_data) +
			(1+1 + 1+1 + 1+1 + 1+1)*blocksize;
	
			
			printf("filesize=%g [MB]\n", filesize/1024.0/1024.0);
			printf("ntot=%u\n", n_tot);
			printf("sizeof(particle_data)=%u\n", sizeof(particle_data));
			
			
			if (n_totvm) filesize+=2*blocksize;
			if (header1.npart[0]) filesize+=4*blocksize;
			if (header1.flag_cooling) filesize+=2*blocksize;
			
			if (mem0) delete [] mem0;
			mem0 = mem = new char[filesize];
			
			printf("filesize[%u]: %u\n", ifile, filesize);
			printf("sizeof(io_header_1): %u\n", sizeof(io_header_1));
			
			*(uint*)mem = sizeof(io_header_1);
			mem+=blocksize;
			memcpy(mem, &header1, sizeof(io_header_1));
			mem+=sizeof(io_header_1);
			*(uint*)mem = sizeof(io_header_1);
			mem+=blocksize;
			
			
			// Pos
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
			for (uint k=0; k<pTypeNum; k++)
			for (uint j=0; j<header1.npart[k];j++) {
				memcpy(mem, P[ this->ofs[k] + ofs[k] + j ].Pos, 3*sizeof(float));  
				mem+=3*sizeof(float); 
			}
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
		
			// Vel
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
			for (uint k=0; k<pTypeNum; k++)
			for (uint j=0; j<header1.npart[k];j++) {
				memcpy(mem, P[ this->ofs[k] + ofs[k] + j ].Vel, 3*sizeof(float));  
				mem+=3*sizeof(float); 
			}
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
			
			// Ids
			*(uint*)mem = n_tot*sizeof(float);
			mem+=blocksize;
			for (uint k=0; k<pTypeNum; k++) {
				memcpy(mem, &Id[this->ofs[k]+ofs[k]], header1.npart[k]*sizeof(uint));  
				mem+=header1.npart[k]*sizeof(uint);
			}
			*(uint*)mem = n_tot*sizeof(float);
			mem+=blocksize;
			
			// Var Masses
			if (n_totvm) {
				*(uint*)mem = n_totvm*sizeof(float);
				mem+=blocksize;
				for (uint k=0; k<pTypeNum; k++) 
				for (uint i=0; i<header1.npart[k]; i++)
				if (header1.mass[k]==0) { 
					memcpy(mem, &P[this->ofs[k]+ofs[k]+i].Mass, sizeof(float));  
					mem+=sizeof(float); 
				}
				*(uint*)mem = n_totvm*sizeof(float);
				mem+=blocksize;
			}
			
			// Gas
			if (header1.npart[0]) {
				
				
				// Internal Energy
				*(uint*)mem = header1.npart[0]*sizeof(float);
				mem+=blocksize;
				for (uint j=0; j<header1.npart[0]; j++) {
				memcpy(mem, &P[*ofs+j].U, sizeof(float));  
				mem+=sizeof(float);	}
				*(uint*)mem = header1.npart[0]*sizeof(float);
				mem+=blocksize;
				
				
				// Density
				*(uint*)mem = header1.npart[0]*sizeof(float);
				mem+=blocksize;
				for (uint j=0; j<header1.npart[0]; j++) {
				memcpy(mem, &P[*ofs+j].Rho, sizeof(float));  
				mem+=sizeof(float);	}
				*(uint*)mem = header1.npart[0]*sizeof(float);
				mem+=blocksize;
				
				// Cooling
				if(header1.flag_cooling) { // this is for the cooling in the user guid it says it is currently unused
					*(uint*)mem = header1.npart[0]*sizeof(float);
					mem+=blocksize;
					for (uint j=0; j<header1.npart[0]; j++) {
					memcpy(mem, &P[*ofs+j].Ne, sizeof(float));
					mem+=sizeof(float);	}
					*(uint*)mem = header1.npart[0]*sizeof(float);
					mem+=blocksize;
				} 
				
			}
			
			printf("file[%u]:\n", ifile);
			char ifilename[256];
			header1.num_files>1 ?
			sprintf(ifilename, "%s.%u", filename, ifile) :
			sprintf(ifilename, "%s", filename);
			
			
			FILE *pFile = fopen (ifilename , "wb" );
			fwrite(mem0, sizeof(char), filesize,  pFile);
			fclose(pFile);
			
			
			for (uint k=0; k<pTypeNum; k++) 
			ofs[k]+=header1.npart[k];
			
			printf("file[%u] is passed.\n\n", ifile);
		}
		
		for (uint k=0; k<pTypeNum; k++)
		printf("ofs[%u]: %u\n", k, ofs[k]);
		mem = mem0;
		delete [] ofs;
	}
	
	
	
	
	void print_P(uint nprint){
		
		for (uint k=0, ofs=0; k<pTypeNum; ofs+=header1.npartTotal[k], k++) {
			printf("%u: \033[22;37m%u, %u\033[22;30m\n", k, ofs, header1.npartTotal[k]);
			
			if (header1.npartTotal[k]) 
				for (uint j=0; j<=nprint; j++){
					uint ik = (header1.npartTotal[k]-1)*double(j)/nprint; // becouse double precision problem
					//uint ik = j;
					uint i = ofs+ik;
					
					printf("\033[22;34mP\033[22;30m[%7u]:  \033[22;35m%10u\033[22;30m  % .7e", i, Id[i], P[i].Mass);
					printf("\t% .7e % .7e % .7e", 	P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
					printf("\t% .7e % .7e % .7e", 	P[i].Vel[0], P[i].Vel[1], P[i].Vel[2]);
					if (k==0 && header1.npartTotal[k]>0) 
					printf("\t% .7e % .7e % .7e", 	P[i].U, P[i].Rho, P[i].Ne);
					printf("\t\033[22;36m %u \033[22;30m", P[i].Type);
					printf("\n");
					
				}
		}
		
	}
	
	void print(){
		
		const double MB = 1<<20;
		
		printf("\n");
		printf("\033[22;34mFilename\033[22;30m: \033[22;32m'%s'\033[22;30m\n", filename);
		
		printf("\033[22;34mFilesize\033[22;30m: %.3f MB \033[22;37m( %u )\033[22;30m\n", filesize/MB, filesize);
		printf("\033[22;34mReed Speed\033[22;30m: %.2f MB/sec\n", readspeed/MB);
		printf("\033[22;34mReed Time\033[22;30m: %.2f sec\n", readtime);
		printf("\033[22;34mNumFilesPerSS\033[22;30m: %u\n\n", header1.num_files);
		
		printf("\033[22;34mBoxSize\033[22;30m: %.15g\n", header1.BoxSize);
		printf("\033[22;34mOmega0\033[22;30m: %.15g\n", header1.Omega0);
		printf("\033[22;34mOmegaLambda\033[22;30m: %.15g\n", header1.OmegaLambda);
		printf("\033[22;34mHubbleParam\033[22;30m: %.15g\n", header1.HubbleParam);
		printf("\033[22;34mRedshift\033[22;30m: %.15g\n", header1.redshift);
		printf("\033[22;34mTime\033[22;30m: %.15g\n\n", header1.time);
		
		printf("\033[22;34mNumPart\033[22;30m: %u\n", NumPart);
		printf("\033[22;34mNgas\033[22;30m: %u\n", Ngas);
		printf("\033[22;34mN\033[22;30m: \033[22;35m%u\033[22;30m\n\n", Ngas+NumPart);
		
		for (uint i=0; i<pTypeNum; i++)
			header1.npart[i] ?
			printf("\033[22;34mnpart\033[22;30m[%u] = \033[22;33m%u\033[22;30m\n", i, header1.npart[i]) :
			printf("\033[22;37mnpart[%u] = %u\033[22;30m\n", i, header1.npart[i]);
			printf("\n");
		
		for (uint i=0; i<pTypeNum; i++)
			header1.mass[i] ?
			printf("\033[22;34mmass\033[22;30m[%u] = %.15g\n", i, header1.mass[i]) :
			printf("\033[22;37mmass[%u] = %g\033[22;30m\n", i, header1.mass[i]);
			printf("\n");
		
		for (uint i=0; i<pTypeNum; i++)
			header1.npartTotal[i] ?
 			printf("\033[22;34mnpartTotal\033[22;30m[%u] = \033[22;33m%u\033[22;30m\n", i, header1.npartTotal[i]) :
			printf("\033[22;37mnpartTotal[%u] = %u\033[22;30m\n", i, header1.npartTotal[i]);
			printf("\n");
		
	}
	
	
	
	
};

#endif


// MODIFICATIONS

/* 1. In version .02 was memory consumption problem after deconstruction of GSIO object.
 * The reason was after deletion on pointers to arrays erased and that is way you cant 
 * to free memeory from prevous using like this: if (p) delete [] p
 * 
 * 2. there was also some problem with built-in baseName Finder, so i replased it with FND,
*/


  











